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Abstract 

The periodic standing wave approach to binary inspiral assumes rigid rotation of gravitational 
fields and hence helically symmetric solutions. To exploit the symmetry, numerical computations 
must solve for "helical scalars," fields that are functions only of corotating coordinates, the labels on 
the helical Killing trajectories. Here we present the formalism for describing linearized general rela- 
■ tivity in terms of helical scalars and we present solutions to the mixed partial differential equations 

D ' of the linearized gravity problem (and to a toy nonlinear problem) using the adapted coordinates 

pL^ | and numerical techniques previously developed for scalar periodic standing wave computations. 

We argue that the formalism developed may suffice for periodic standing wave computations for 
post-Minkowskian computations and for full general relativity. 
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The computational study of the inspiral of binary black holes is important for the understanding of gravitational 
wave signals, and is of inherent interest as a question in general relativity that can be answered only with computation. 
It has therefore become the focus of supercomputer codes that evolve Einstein's field equations forward in time 
from initial conditions chosen to represent a starting configuration of the inspiralling objects. The evolution codes, 
however, typically become unstable on a timescalefset by the size of the hole) short compared to a full orbit. Reliable 
calculations of the final plunge are now feasible the merger and ringdown of the final black hole fate of the 
system are handled well with perturbation theory Q, and the early inspiral is well approximated with post-Newtonian 
£jTj! computations^]. What cannot be handled well is the intermediate phase of the inspiral, the late epoch during which 
• • nonlinear effects are too strong for a post- Newtonian approximation, but for which too many orbits remain for stable 
. £h ' numerical evolution. 

It has long been recognized that the basis of an approximation scheme should be the slow rate of inspiral, the small 
'_ ■ ratio of the orbital time to the radiation damping time0,[(|. Through an adiabatic treatment of the slow inspiral, such 
an approximation could give answers about the radiation and rate of inspiral in the intermediate epoch. In addition, 
when the rate of inspiral becomes too rapid, the intermediate approximation could hand the problem off to numerical 
evolution codes to do the final orbit or two and the plunge. Along with the problem being handed off, would be 
the ideal initial data for the subsequent evolution. The need for and the concept of an intermediate approximation 
have been clear, but such an approximation has not been easy to implement. Along with several others0> H, 0, 
we have based an approximation for slow inspiral on a numerical computation of no inspiral. That is, we seek a 
numerical solution of Einstein's equations for binary objects that are in circular periodic motion, and whose "helically 
symmetric" fields rotate rigidly with the source objects. (For a definition of helical symmetry, see Sec. II of Ref. [Tl|.') 
The universality of gravitation suggests that the unchanging motion of such a system is not compatible with 
outgoing radiation, and this intuitive suggestion is confirmed by the mathematics of the theory. We therefore seek 
a helically symmetric solution for the sources coupled to standing waves, not to outgoing waves. In a linear theory, 
standing waves, in the sense that we use the term, are a superposition of half-ingoing and half-outgoing solutions. 
From the fact that the solution, in linear theory, is half the superposition of the ingoing and outgoing solutions, one 
could extract the outgoing solution. The crux of our periodic standing wave (PSW) method is that even for highly 
nonlinear binary inspiral fields there is an "effective linearity." The standing wave solution, to good accuracy, is half 
the sum of the outgoing plus ingoing solutions despite the nonlinearities. In general relativity, therefore, our goal is 
to solve the standing wave numerical problem and from that solution to extract an approximation to the outgoing 
solution. 

It is important to understand why effective linearity can be correct for inspiral. In the strong-field regions very 
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close to the sources, the solution is very insensitive to the choice of distant radiative boundary conditions (ingoing, 
outgoing, standing wave) . In this near-source region a superposition of half the ingoing and half the outgoing solution 
gives a good approximation solution, because it amounts to averaging two samples of the same thing. In the wave 
zone where the outgoing and the ingoing solutions are very different, the fields are weak enough that nonlinear effects 
are negligible, and once again we can superpose. The strong-field region and the boundary-influenced region should 
be widely separated unless the sources are rotating very close to c, in which case the wave zone will start just outside 
the sources. It is, however, not expected that ultrarelativistic source motion can occur during the slow inspiral epoch 
of motion. 

We have recently ^3 been able to confirm effective linearity. This confirmation has been achieved with a model 
problem, since the validity of effective linearity can only be carried out in a model problem. In general relativity, there 
can be no "true outgoing" solution available for confirmation until numerical evolution codes are fully developed. In 
addition, the numerical features of the helically symmetric standing wave calculation pose new challenges very different 
from those of evolution codes, and are best resolved in the simplest context possible. 

The model problem in Ref. 0] used a nonlinear scalar field theory with a pair of diametrically opposite point-like 
scalar charges. The imposition of helical symmetry on the problem leads to a boundary value problem for a system 
of mixed (hyperbolic and elliptical) partial differential equations. To solve that boundary value problem efficiently, a 
set of numerical techniques was developed that we called the "eigenspectral" method. In the present paper we report 
two important steps toward using the PSW method for full general relativity: First, we develop the infrastructure 
for describing gravitational fields. In previous work with scalar toy models helical symmetry was imposed by simply 
requiring that the scalar field is a function only of three corotating coordinates (labels on the Killing trajectories). 
We find it useful to use the expression "helical scalars" for such functions that depend only on corotating coordinates. 
Our computation is done on a grid of corotating coordinates, so straightforward computations can only be carried out 
for helical scalars. But complications arise with tensor fields. The components of a helically symmetric tensor field 
are generally not helical scalars. A resolution of this difficulty is to compute only with projections of the tensor on 
a helically symmetric basis, that is, on a basis that is Lie dragged by the helical symmetry. These projections would 
be helical scalars. A helically symmetric basis, however, has its own dynamics, and complicates the time dependence 
of the projected fields. In the infrastructure developed here we show that the use of a corotating "pure-spin" basis 
leads to a remarkably simple description of helically symmetric tensorial fields. That infrastructure is presented in 
the explicit context of linearized general relativity. 

The second step taken in this paper is to present numerical results showing that no new computational problems 
are encountered in dealing with the helical scalars of the linearized gravity problem. Since no new problems were 
anticipated this numerical work simply constitutes a confirmation that the "eigenspectral method," the set of tech- 
niques developed for scalar models, appear to work equally well for linearized gravity. Those techniques include the 
use of (i) "adapted coordinates," a corotating coordinate system that conforms well to the source regions and to the 
radiation field, (ii) "multipole filtering," the elimination of numerical noise associated with angular differencing, by 
keeping only a few multipoles of the adapted coordinates, and (iii) the modification of the multipoles so that they are 
computationally orthogonal at the level of machine precision. In this paper we provide numerical solutions using this 
set of techniques. Since these solutions differ very little from the numerical problems studied in detail in Ref. [T^ 
the presentation of linearized results is brief. The solutions to nonlinear problems are much more difficult than those 
for linearized problems. Newton-Raphson iteration must be used, and convergence of the iterative process has been 
the major challenge in numerical work. Again, there is no apparent reason the problem should be more difficult for 
the tensor-based helical scalars than for the nonlinear scalar models of Ref. [l^- Nevertheless, it is an issue worth 
checking, and initial results are briefly reported showing that the scalar techniques work well for a model tensorial 
problem with a simple toy nonlinearity. 

The rest of this paper is organized as follows. Section ^ gives a brief review of the scalar PSW problem in order 
to introduce adapted coordinates and the application of multipole filtering and the eigenspectral method in those 
coordinates. Section ITTll presents the description of the helically symmetric fields of linearized gravity that is suitable 
for computation. The field equations are given for general corotating coordinates and series solutions are given, in 
corotating spherical coordinates, for the problem of two equal masses in circular orbits around each other. The field 
equations for linearized gravity using adapted coordinates are given in Scc. lIVI In that section, also, are given the forms 
in adapted coordinates of other elements of the computational problem, the inner boundary condtions specifying the 
sources, and the outer boundary conditions specifying the nature of radiation. Section presents numerical results, 
comparing the series solutions of Sec. IIIII with the solutions of the field equations using the eigenspectral method. In 
addition, in this section a description is given of a toy nonlinear modification to linearized gravity, and results are 
given demonstrating that the resulting "theory" can be solved with the techniques that worked for the scalar case. 
The implications for the next steps in the PSW program are discussed in Sec. I VII In particular, it is argued that 
almost all the infrastructure for solving the PSW problems in the post-Minkowski approximation and in full general 
relativity problem is probably established in the work with the linearized problem. The Appendix gives some detailed 
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expressions needed for computations in adapted coordinates. 

Throughout the paper we use units in which c = G = 1 and we follow the conventions of Misner, Thorne and 
Wheeler 



II. SCALAR MODELS, COORDINATES, AND NUMERICS 

A. Nonlinear scalar models 

Here we review the model problem and numerical techniques of Ref. |T^ | . Our model problem was a nonlinear scalar 
field coupled to point-like sources in Minkowski space, and satisfying the held equation 

* ;a; /33 a/3 + F = V 2 * - -^8^ + F = Source . (1) 

In principle, the source was taken to be two points of unit scalar charge in orbit around each other at angular frequency 
O, and at radius a. In practice, the computational problem used inner boundary conditions on small, approximately 
spherical surfaces to represent the effect of a point source; no source term was included in the held equation that was 
computationally solved. The velocity parameter for the system v = aft was taken to be of order unity, representing 
the strong-held tight binary for which post-Newtonian approximations are inadequate. 

The term F contains the nonlinearity in our scalar model theory, and we found the following form, with parameters 
A and to be useful: 

A * 5 

A crucial feature of F is that, like the nonlinearities of general relativity, it is very large near the sources, and becomes 
negligible far from the sources. The A multiplier allowed us to vary the strength of the nonlinear term, and the 
parameter allowed us to vary the profile of the nonlinearity in the strong held region. 

Our scalar problem was defined by Eqs. and (0), and by the source motion at angular frequency ft in the 
equatorial plane. With standard spherical coordinates, helical symmetry can be imposed on the solution ty(t,r,9,(j)) 
by restricting to solutions of the form ^(r, 9, ip), where (p is the comoving azimuthal coordinate cj>— ftt. By restricting 
the solution in this way, we have eliminated the possibility of "evolution." For such helically symmetric solutions a 
change in time by At is the same as a change the azimuthal angle A(j> = — ftAt. 

When the restriction to helical symmetry is made, the held equation becomes 

= — — (r 2 — \ + 1 9 ( S [ n 0^1 
r 2 dr V dr ) r 2 sin 9 88 \ 89 



r 2 sin 2 i 



ft' 



8 2 V 

dtp 2 



Source - = ct(#) , (3) 



and the mixed nature of the partial differential equation becomes obvious. The principal part of this quasilincar 
equation is elliptic inside a cylinder at rsin# = l/ft, and hyperbolic outside that cylinder. The problem is to be 
solved with radiative conditions (ingoing, outgoing, or standing wave as described below) on a spherical surface at 
large distances from the sources. Well posed problems in physics typically supply cauchy data on open surfaces to 
hyperbolic equations, and Dirichlet or Neumann data on closed surfaces to elliptic equations. Our model is unusual 
in that it leads to a boundary value problem with "radiation" conditions on a closed surface surrounding a mixed 
problem. Though unusual, our problem is intuitively well posed, and passes a computational test: we have found 
no fundamental difficulty in solving models of this type numerically. Furthermore, a careful analysis^) proves that 
solutions exist and are stable for a closely related problem. 

"Standing wave" solutions - half ingoing and half outgoing - are at the heart of our method, but there is not an 
unambiguous definition of standing wave solutions in a nonlinear theory. Our procedure is to hnd the outgoing C~^ t 
and ingoing C^l Green functions for Eq. @. In principle, we can then iterate to hnd a solution of Eq. @. The 
iteration 

, (4) 
if it converges, gives ^out, our nonlinear outgoing solution (and similarly for ^m), while the convergent result of 

^^Ite+A^^tl)] (5) 
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is what we mean by our nonlinear standing wave solution 'J'std- The standing wave solution 'J'std is fundamentally 
different from (\I/ ou t + ^ , i n )/2, but if effective linearity is correct, the two are very nearly equal. (Note: In practice, 
for strong nonlinearities, the direct iteration described above must be replaced by Newton- Raphson iteration.) 

A central idea of the PSW approximation is that the "exact" (i.e. , numerical) solution to the standing wave 
problem is an excellent approximation to half the sum of the outgoing and ingoing solutions. If this is so, it means 
that from the standing wave solution we should be able to extract an excellent approximation to the outgoing and the 
ingoing solutions. (It should be noted that this statement is meaningful for nonlinear model field theories in which 
there is meaning to an outgoing helically symmetric or ingoing helically symmetric solution to the field theory. For 
general relativity, a helically symmetric spacetime with outgoing radiation is impossible. As explained in Ref. |12| . 
the relevance and justification of the PSW approximation for full general relativity lie in the fact that the method 
gives an approximation only for a limited region of spacetime.) 

The extraction of the approximate nonlinear outgoing solution from the computed standing wave solution is a 
direct application of the concepts underlying the argument for effective linearity. In the weak wave zone, far from the 
sources, the field theory is very nearly linear and the solution must be very nearly a standing wave solution to that 
linear theory. It is, therefore, straightforward to deconstruct it into outgoing and ingoing solutions. The extracted 
outgoing solution can be continued inward through the induction zone into the near-field zone. If the theory were 
completely linear, this continuation would be valid up to the source points. Due to the nonlinearity, however, this 
procedure is no longer valid at distances so close to the source points that the nonlinearity of the theory is important. 

In the near field region close to the source points, at a distance small compared to a wavelength, the solution is 
highly insensitive to the nature of the boundary conditions (i.e. , whether they are outgoing or ingoing). Here, we can 
use the standing wave solution itself as the outgoing solution. This near field solution should extend out to the weak 
field region and overlap with the extracted outgoing solution described above. These two solutions, the weak field 
outgoing solution outside the strong field region, and the standing wave solution in the strong field region, are then 
patched together, with some blending in the region of overlap, and the result is the PSW extracted approximation to 
the nonlinear outgoing solution. 

B. Adapted coordinates 



z Y 




x Z 

FIG. 1: Two systems of corotating Cartesian coordinates. 

It is useful to identify a number of coordinate systems to describe the physical problem, including three Cartesian 
systems. The system {x, y, z} is that of inertial Cartesian system in which the z axis is the rotation axis for the source 
objects, with rotation in the positive sense about the z axis. In general we use tildes to distinguish the corotating 
version of a quantity when that quantity occurs in forms both related to inertial and to corotating systems. The 
set {x, y, z}, in this sense, is the corotating set of Cartesian coordinates for which z = z and for which the source 
points remain fixed on the x axis. The system {r, 9, (/>} is that of inertial spherical coordinates defined in the usual 
way relative to {x,y,z}. The system {r, 9, ip} is a set of coroting spherical polar coordinates, defined by the usual 
transformation relative to {x, y, z} . The two systems of spherical coordinates are related by (p = <j) ~ The 
Cartesian system {X, Y, Z} is a convenient renaming of {x, y, z}, with Y the rotational axis, and Z 7 the axis through 
the source points. 

Our adapted coordinate system {\, 6, *&} is a corotating two-center bipolar coordinate system defined, relative to 
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the {X, Y, Z} system, by 



(z-a) 



X 2 + Y 2 



Z + a) +X 2 + Y 2 



e 



1 / 2Z\/X 2 +Y 2 

- tan — — 

2 \Z 2 -a 2 -X 2 -Y\ 

tan" 1 fy/x) , 



1/4 



(6) 

(7) 
(8) 



where a is the distance from the center to each of the source points. 

For x 3> °i the adapted coordinates become a corotating system of spherical coordinates defined relative to 
{X,Y,Z}. (That is, the Z axis is the azimuthal axis for the $ coordinate.) The adapted coordinates, pictured 
in Fig. |3 are discussed in greater detail in Ref. |T^. 




FIG. 2: Two-center bipolar adapted coordinates. On the left is shown curves of coordinates \ an d in the $ = orbital 
plane. The surfaces on the right show surfaces of constant x> ©, and 3>. Here a is the distance from either of the two centers 
(locations of the point sources) to the midpoint between the centers. The x, y, z coordinate system is a corotating Cartesian 
system for which the z axis is the rotation axis; for the X, Y, Z a corotating Cartesian system the Z axis is the line through 
the centers. 

We originallvfToj solved the model problem of Eq. © in the {r, 9, tp} system with more-or-less straightforward finite 
differencing and direct matrix inversion. (The mixed nature of the partial differential equations prevents the use of 
such efficient techniques as overrelaxaton.) This approach was successful (iterations converged) for models with a 
limited range of source velocities and nonlinearities. Subsequently we developed an innovative numerical method that 
gave remarkably good results for the scalar problem, with very little computational cost. The new method is based 
on three elements. 

First, we used the adapted coordinates discussed above. The solutions to the finite difference form of Eq. in 
these coordinates turned out to be plagued by what appeared to be angular noise. This noise was eliminated with 
multipole filering, a form of smoothing of the angular variation. In this method the scalar field was expanded in 
spherical harmonics of the angular functions and <&. The fact that the adapted coordinates are well suited both to 
the source structure and to the radiation field suggests that good accuracy can be achieved when only a few multipolcs 
are kept. In fact, good accuracy was found when only the monopole and quadrupole were kept in the case of sources 
speeds around 0.3c -0.4c or less. At higher speeds the radiation field develops sharper gradients and more multipoles 
must be kept to achieve reasonable accuracy. 



C. inner and outer boundary conditions 



For identical point sources of unit scalar charge moving in the equatorial plane at angular velocity Q in circular 
orbits of radius a, the source used in Eq. Q was 

S = - t- 1 ^^ 6(6 - tt/2) [%) + <% - tt)] , (9) 



with 7 = 1/ VI — v 2 =1/ vl — a 2 fl 2 . For this source, it was shown in Ref. [l2j that the small- \ limit of $ is 

1 2a 1 
4ttx 2 y/i + 7 2 w 2 S i n 2 20 cos 2 $ 



* = -~^= . ; . ,^ ^ ■ (10) 
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We take this as an inner boundary condition at a small value of \ that determines one limit of our computational 
grid. 

The outer boundary condition used is simply the ingoing or outgoing Sommerfeld condition expressed in terms of 
the adapted coordinates. In Ref. 01 it was shown that in adapted coordinates this condition becomes 

in which the upper and lower signs correspond, respectively to the outgoing and ingoing conditions. 



D. multipole filtering and the eigenspectral technique 

In practice, multipole filtering was carried out as follows. If there are N grid locations @i, <F,, of the angular 
coordinates 0, $, then at each value of the radial coordinate \ there are N values of the scalar field ^f(x, ©i, ln 
multipole filtering, a set of M < N spherical harmonics is used, at each value of x> as weighting factors for these 
scalar field values, and M weighted sums (i.e. , multipole projections) are taken of the scalar field values. In the same 
manner, the N field equations at x are projected into M sums. Those M projected equations arc then solved for the 
M projections of the field. 

It was found that this procedure did not work if the multipole weights Yy were simply taken to be YLm(Oz, ^j), 
the continuum spherical harmonics evaluated at the discrete angular grid locations. These projection weights are only 
orthogonal in the continuum limit, and their failure to be numerically orthogonal to high precision was the probable 
source of angular noise that plagued computations. Slightly modified weights, in place of the grid-evaluated spherical 
harmonics, gave us weights that was orthogonal to the level of machine precision. The use of these modified weights in 
multipole filtering eliminated the problem of angular noise, while at the same time significantly reducing the number 
of equations to be solved, and hence reducing the computational burden. See Ref. |12j for details. 

We focused, in Ref. [12j, on the most important question that can be answered with these models and numerical 
methods: Does effective linearity work? Can we extract a good approximation to the outgoing nonlinear problem 
from the sort of standing wave computation we will be limited to when dealing with Einstein's theory? Figure |3| 
gives strong evidence that we can. In that figure, the computed outgoing nonlinear solution is shown as a solid curve. 
The data-type points represent the outgoing solution extracted in the manner described above in Sec. Ill Al For the 
parameters A = —15 and *S? = 0.15 in Fig. El nonlinearities are significant, strong enough to reduce field strength by 
around two-thirds. The outgoing and standing wave solutions were each computed by the Newton-Raphson version 
of the iteration in Eqs. Q),j5j|. We have run models with much stronger nonlinearity and have found equally good 
agreement of the true outgoing solution and the extracted approximation. The validity of effective linearity should, 
in fact, become questionable not for stronger nonlinearity, but only for physically implausible high source velocity. 

In addition to confirming effective linearity, computation with the model has also allowed some early insights about 
sensitivity to source details. By varying the multipole content of the inner boundary data we explored the impact of 
source structure on the radiation field. The result (detailed in Ref. 0|) is in perfect accord with physical intuition; the 
radiation is insensitive to source structure unless the source size becomes comparable to the separation of the sources 
(i.e. , unless the moments ascribable to the structure of the individual sources are comparable to the quadrupole 
moment due to the separation of the mass points). The equivalent question for Einstein's theory is more difficult, but 
we should eventually be able to give clear quantitative answers. 



III. THE DESCRIPTION OF HELICALLY SYMMETRIC LINEARIZED GRAVITATIONAL FIELDS 



A. The physical problem 



We use the standard description and notation for the linearization of Einstein's field equations (see Chap. 8 of 
Ref. jl3j|). The perturbations from flat spacetime /i^„ are defined by 

9nv = "r\^v + h^ v , (12) 
and the trace-reversed perturbations h^ v are defined by 

hfj, u = hfj, u - -rjfj, v h , (13) 



FIG. 3: The computed nonlinear outgoing solution compared with an approximate outgoing solution extracted from the 
computed nonlinear standing wave solutions. The solid curves show the computed outgoing solution for a scalar model with 
nonlinearity parameters A = — 15 and — —0.15 and with source velocity aQ — 0.4c. The data-like markers show the 
approximate outgoing solution extracted from the standing wave solution; the shape of the markers indicate whether the 
extracted solution corresponds to the deconstruction of the linearized solution, the standing wave solution, or a region of 
blending of the linearized outgoing and standing wave solutions. The computations of both the standing wave solution and the 
outgoing solution were carried out using quadrant symmetry with a 40 x 80 angular grid, and with 16001 divisions in \ ranging 
from x — 0.02a to 80a. The multipole filtering kept only the multipoles corresponding to £ = 0, 2 and 4. 

where h = rf^h^y. For linearized computations, is treated as a tensorial field in Minkowski spacetime, and indices 
are raised and lowered with the Minkowskian metric rj^ v . In the gravitational Lorentz gauge 

h% = 0, (14) 

the linearized field equations of general relativity are 

V,a Q = -16^. (15) 

Our physical problem is that of two points, each of mass mo, in circular orbits , with radius a, angular velocity f2, 
and hence speed v = aQ through the background Minkowski spacetime. The general form of the stress-energy source 
for linearized theory [l5| is 



J—oo CLT ar u r 



in which 



(fl-fi(i)), (16) 



U u = 7 = l/Vl-v 2 U x = TwysmSlt U y = ±vjcosnt. (17) 



The signs of U x and U y are different for the two source particles. One, call it particle 1, is at <f> = £lt; the other, 
particle 2, is at cf) — fit + ir. In Eq. 117(1 the upper sign corresponds to particle 1, the lower to particle 2. The explicit 
nonvanishing components, in the inertial t, r, 9, <j) system, are then, 

T tt = m 7 S{r ~ a) S (9 - tt/2) [S (if) + 6 {<p - tt)] (18) 

T tx = TvsinntT u T ty = ±v cos Clt T u (19) 

T xx = v 2 sin 2 fit T u T yy = v 2 cos 2 fit T u T xy = -v 2 sin fit cos Qt T u . (20) 

We use this helically symmetric stress energy as an explicit source when we derive the series solutions to Eq. <|15H • 
For the solution of Eq. I|15|) via the eigenpectral method of Ref. |l2|. however, we find hap in the limit of small distance 
from the particles, and impose this solution as inner boundary conditions on the homogeneous form of Eq. (|15|l . To 
get the near-particle solution we use the Lienard-Wiechert solution of Eq. (|15[1 



(21) 
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The retardation condition "ret," means that r is to be evaluated at time such that 

\x-f(r)\ = 0. 



(22) 



For the particle at (p = 0, we now evaluate this approximately by assuming that the particle moves with constant 
velocity x = a, y = aClt pax t = i>i par t, z = 0, through the inertial coordinate frame. At field point t, x, y, z we have the 
following retardation condition for the particle at x — a: 



= -(t - ipart) 2 + {x- a) 2 + (y- vt palt ) 2 + z 2 , 



with the solution 



Cpart — 7 



t — yv — y/ (t — vy) 2 7 (r 2 — t 2 )/j 2 
in which 

2 / \2 , 2 , 2 

r — (x — a) + z + y . 
To express this in corotating coordinates, we next introduce the notations and approximations 



~ . . . - - ~2 /— \2 , ~2 , ~2 

y = y + vt x = x z = z r = (x — a) + y + z 



With these, Eq. QHty simplifies to 



kpart 



t — v^f 2 y — ^yjr 2 7 7 2 w 2 y2 , 



(23) 



(24) 



(25) 



(26) 



(27) 



and, finally, we can evaluate 



U ■ (x-r{r)) =-U°(t- < part ) 7 U y {y - vt paxt ) 

ret 

= -j[t- t part -v(y- vt pait )} 



= - 7 1 (t ~ *part) + "72/ = - Vr 2 + J 2 v*y 



2 712 



(28) 



With Eqs. HI 7(1 and ((28(1 . the expression in Eq. ((21(1 gives us the inner boundary conditions to be used for solving 
Eq. JTSJ: 



= 4m 

h nx = h xn = 4m 

/i na = /i yn = 4m 

= 4m 

h w = 4m 

^xy = - h yx = 4mQ 





+ 7 2 v 2 y* 2 
7 2 




_|_ r^-ylyl 

7 2 




7 7 2 u 2 j/ 2 
7 2 




_|_ ry2 v 2yl 

7 2 


x/r 2 


7 r ) 2 v 2 y 2 
7 2 


V? 2 


7 7 2 u 2 j7 2 



[=Fu sin fit] 

[±i> cos fit] 

[u 2 sin 2 fit] 

[w 2 cos 2 fit] 

[— f 2 cos fit sin fit] 



(29) 
(30) 

(31) 
(32) 

(33) 
(34) 



The outer boundary conditions, roughly speaking, are the conditions that nonradiative parts of the field fall off as 
1/r™, and the radiative parts of the field satisfy simple ingoing or outgoing Sommerfeld conditions. The details are 
given after the presentation of the formalism for describing helically symmetric tensorial fields. 
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B. Description of helically symmetric tensorial fields 

Imposing helical symmetry on the tensorial field means that the Lie derivative of h^ v vanishes along the helical 
Killing field, or 

MV) = o. (35) 

Here the helical Killing vector is 

£ = dt + nd4 = dt. (36) 

The first expression gives £ in terms of the inertial spherical background coordinates t, r, 6, </>, and the second in terms 
of the corotating spherical background coordinates t,r,6,Lp — tf> — Qt. As explained in Sec. [IJ our computational 
unknowns are fields on a grid of corotating coordinates. We must therefore cast the field equations in terms of a set 
of functions that are "helical scalars," functions only of the spatial corotating Minkowski coordinates such as {r, 0, (p} 
{x, y, x}, {X, Y, Z}, or the adapted coordinates {%, 0, $}. 

The key to this is to project on a corotating basis. A corotating vector V is one for which £g(V) = 0. If V and W 
are both corotating vectors, and is helically symmetric, then h^V^W 1 ' is a helically symmetric scalar, a function 
only of three corotating coordinates. Our approach, then is to use a coroating basis and to project all components of 
h^ v on this basis. 

If not done with some care, projection on a corotating basis, can destroy the simplicity of the linearized field 
equations (|15|l . Since the inertial basis is covariantly constant, the equations separate for the components of h^ v ; each 
component satisfies its own equation, and the system of equations separates into a set of scalar-like equations. This 
is not true in general for the scalars formed by projections with the corotating bases. Choices can be made, however, 
that result in a high degree of separation, and a very simple set of equations. 

To achieve this simplicity we start by defining the covariantly constant orthonormal inertial basis system for the 
Minkowski background: 

n = d t e x = d x e y = d y e z = d z . (37) 

Here we closely follow the analysis given by Thorne^fJ. That analysis uses the method of Mathews |l7|. and introduces 
a set of second rank symmetric spatial basis tensors t2j of spin 2, i.e. , tensors that transform among themselves as 
an irreducible represenation of the rotation group of order 2. To these we add an additional spin spatial tensor to,o 
and the spin and spin 1 second-rank symmetric tensors needed to include the timelike direction 

(38) 
(39) 
(40) 
(41) 
(42) 
(43) 
(44) 

We next define the corotating equivalents to the basis vectors 

n = n e x = e x cos fit + e y sin ilt e x = —e x sin fit + e y cos fit e z = e z , (45) 
and we use these to define the corotating equivalents of the basis tensors: 

(46) 
(47) 
(48) 
(49) 
(50) 
(51) 
(52) 



^nn — 


nn 






t n = 


ne z 






tn,±l = 


f [n(e x 


■ ie y ) 4 


- (e x ± ie y )n] 


to,o = 




+ e y e y 


+ e z e z ] 


t2,0 = 




+ e y e y 


- 2e 2 e z ] 


*2,±1 = 


~^~2 [^x^z 


+ e z e x ] 


2 ^ ~t~ 


t2,±2 = 


2 


- e y e y ± 


% (& y Q x ~\- G X G. 



t»nn — 


n I' = t nn 




tnO = 


ne z = t„o 




tn,±l = 


±± [n(e x ± ie y ) + (e x ± ie y )n] = 


e* ia H n<±1 


to,o = 


-7= [e x e x + e y e y + e z e z ] = t ,o 




t2,0 = 


y= [e x e x + e y e y - 2e z e z ] = t 2 ,o 




t2,±l = 


~^~2 [^x^z 4~ G^e^] 2 ^ ~\~ ^z 


e y ] = e^ n *t 2l ±i 


t2,±2 


~2 l^x^x ^y^y ^ (.^y^x ^x^y) 


„T2j:nt + 
— e t2,±2 ■ 
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Note that all the tensor bases are orthonormal in the following sense: 

(tij)* ■ (tfc p ) = (tij)* ■ (tk, p ) = SikSjp , (53) 
and that the tensor h can be written as 

h = * ( " n) t ml + * ( " 0) t„ + * (nl) t nl + tf^'-^tn,-! + * (00) t + * (20) t 20 + 

^(21) t2i + ^(2,-l) t2 _ 1 + ^(22) t22 + $ (2,-2) t2 _ 2 f (54) 



or as 



h = ¥ nn h nn + ¥ n0 h n0 + * (lll) t„ 10 + *("'- 1 °)t„,_ 1 + * (00) t o + * (20) t 20 + 

$ (21) t 21 + *( 2 '- 1 )t 2 ,_ 1 + $( 22 >t 22 + $( 2 <- 2 >t 2 ,_ 2 . (55) 

From Eq. (|53l) we get 

■^(nn) = ~ h tt (gg) 

^r(nO) _ ffz ( 5? ) 

^(nl) = _- h tx +i J l ty ^ 

*(°°) = 1 [ft** + + (59) 

*( 2 °) = -=i [h xx + h yv - 2h zz ] (60) 

^(21) _ _J l xz +i J l yz ( 61 ) 

#( 22 ) = I [h** - ftw] _ ihvx _ ( 62 ) 
Here and below, we have dropped functions that are redundant due to the relations 

$K-1) = _ $(2>-l) = - (V 21) )* *(2.-2) _ ^(22) y _ (g 3 ) 

These relations are also true if tildes are placed over all variables. From the rotation behaviors of the t basis, given 
above, the rotation laws for the * fields are: 

(64) 
(65) 
(66) 
(67) 
(68) 
(69) 
(70) 

where and are real functions. 

We define projections of the stress energy by analogy with the projections of the h perturbations, 

yj-(nn) fj-itt ^j-{nn) (71) 

*j~(nO) rptz <j-(nO) (72) 

J-(nl) _ _rptx _|_ ^ rpty _ e -intj-{nl) 

7"(00) _ _J_ !rp X X _|_ rpyy _j_ rp ZZ ^ _ J- (00) (J 

v3~ 

7^(20) — — r^zic _|_ j-yy _ 2T ZZ ] — T^ 20 ^ (75) 

T (2i) _ _ T ^ + irs/z = e -intj-(2i) ( 76 ^ 

7"(22) _ 1 [J 13 ™ _ rpyy^ _ ^ rpy X _ £ - 2 iOt^( 22 ) (jrj^ 



§(«») 


_ ^y(nn) 






$(™0) 


= ^(n0) 






$(«!) 




f/("l) . 


hi y(»i) 


$(°0) 


= #(°°) 






$(20) 


= ^(20) 






$(21) 


= e «M ^(21) = 


C/(21) 4 


. iF (21) 


$(22) 


_ e 2iUt ^(22) _ 


= [/(22) 
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Since the inertial basis tensors are covariantly constant, we can write the field equations (f 1 51) as 

= -WttT^ (78) 
where □ is the simple scalar d'Alembertian. From the relations in Eqs. |(S1J— (J70J and (|71[) -l|77 [l we can then write 

e iknt D f e -iknty(ab)\ = _ 167r <f (ah) . (79) 

with k = 0± 1,±2. 

Lastly, for a helical scalar / we have that 

e lfcm D (e~ lknt f) = □/ + 2iknd t f + k 2 n 2 f . (80) 

The time derivative d t f here uses the time coordinate of the inertial t, r, 9, <fi system and, for / a (helically symmetric) 
function of the corotating coordinates t, r, 6, tp, is equivalent to —Qd v f. The field equations for the metric perturbations 
then take the form 

n¥ ab) - 2ikn 2 o^ + k 2 n 2 ^ = -i67rf (ab) (8i) 

where □ is the scalar d'Alembertian. The explicit equations are 

□$(oo) = _ 16n fm (82) 

□$(n0) = _i 67r r ( ™ 0) (83) 

□$(20) = _i 67r r( 2 °) (84) 

□$(««) = _i67rr ( ™ n) (85) 

□$(«i) _ 2in 2 d ip ¥ nl i + n 2 ¥ n ^ = -Wnf^ (86) 

□$< 21 ) - 2iil 2 d v ¥ 2 V + ft 2 $< 21 ) = -WnfW (87) 
□*( 22 ) - 4iil 2 d v ¥ 22 ^ + 4fi 2 *( 22 > = -167rT( 22 ). (88) 

From Eqs. (|3"4*j) . and the prescriptions in Eqs. ifBfT)) - lfrT2f . and Eqs. f70"|) . we have the inner boundary 

conditions, to be applied at small \ to represent the near field of the mass points at x = > 

= 4m = 72 = ==== (89) 
$(°°) = 4r» ' J_ (90) 



7 2 




_|_ rylylyl 

7 2 




■+- "/ 2 v 2 y 2 
7 2 




_|_ ^2y2y2 

7 2 


^Jr 2 


_|_ ry2y2y2 

7 2 


V? 2 


■+- "/ 2 v 2 y 2 



3 

$(20) = 4 7 = / _ J_ j (Q1) 



V™ = 4m = / = o = t;sgn[cose] (92) 

^(22) = 4m , 7 ( - ^ (93) 

•^/r 2 4- 'y 2 v 2 y 2 \ 2 / 

$(n0) = jj(21) _ y(21) _ ^(nl) = y(22) = q _ /g^ 

C. Outer boundary conditions for helical scalars 

The quantities $( nn ), xpt 00 ), \l>(" ) j $( 20 ) ) are nonradiative multipoles that fall off as 1/r. The radiative parts of 
these fields satisfy a Sommerfcld condition 
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in which the upper and lower signs correspond respectively to outgoing and ingoing waves. In practice, special 
attention is not necessary for the nonradiatable multipoles of ty( nn \ \I>( 00 ) ; \l>(™ ) j i|f( 2 °). The Sommerfeld conditions 
can be applied to the total field, without regard to multipole content. 

Some care must be taken with the Sommerfeld conditions for \j>( nl ), vp( 21 ) and vp( 22 ). The Sommerfeld condition 
applies to the inertial projections of h a p, and hence to ^^ nl \ vI/( 21 ) and \I>( 22 ). With the relations, in Eqs. 1(64(1- JTDJ, 
between these quantities and the "helical scalars" \I>( nl ), \$i( 21 ) : \j/( 22 ) used in computation, we arrive at the conditions 

(96) 



(97) 
(98) 

(99) 
(100) 



d 
dr 




dip 


re i<Pty(nV) 


d 
dr 


re W>vj/( 2 l) 


dip 




— |re 2i¥, $ (22) " 
dr L 


= 

dip L 





or 



in which k = 1 for (a6) = (nl) or (21), and fc = 2 for (ab) = (22). 



r dr V / \ ay 

r dr V / \ oy 



D. Series solutions 



To develop series solutions of Eqs. 1)82(1 - 1)88(1 we start by using the expressions in Eqs. ((18(1 - 120(1 . for stress-energy 
components of the symmetric pair of particles, in the general expressions in Eqs. ((71(1 - ((77(1 . The nonvanishing results 
are 



f(nl) 

7"( 2 0) = U _q-(nn) 



m 7 2 a) S (8 - it/2) [8 (p) + 8 (p - tt)] 



= ±i vT^ 



2 

2 



j<22) = 



(101) 
(102) 
(103) 

(104) 



As in Eq. I(19(l . the upper sign in the expression for jK™ 1 ) indicates the particle at ip = 0, the lower sign indicates the 
particle at (p = tt. 

With these expressions, Eqs. ((8211 ((88|) take the explicit form 



$r(n0) _ q 



□$(00) = -l6^m ^ (r Q2 a) J(cos 0) ^2 + S(ip - tt)} 



(105) 
(106) 
(107) 



□*( 2 °) = + lQTrmo 6 ^ „ a) (5(cos0) [%) + % - *0] 
a V6 



(108) 
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- 2in 2 d v ^ + n 2 ^ = - * i6Trm Q 6{ - r , a) 



<5(cos 9) [S((p) — 6(<p — 7r)] 



(109) 



$(2i) = o 



(110) 



□tf < 22 > - iitfd^W + 4 fi 2 vI/(22) = + Wnm 6ir „ a) ^(cos0) ^ [%) + % - *)] 



(111) 



where 



1 d / 5*\ 1 9 / „<9* 



r 2 Sr \ dr j ' r 2 sin 90 \ " 89 j ' \r 2 sin 2 i 



a 2 * 

<9<p 2 



(112) 



For ifC 20 ), \p( 00 ) the equations all have the form 

_5(r — a) 



= K- 



6 (cos 9) [5(tp)+5((p- tt)] 



(113) 



in which the value of if can be read from Eqs. (|1U5|I - l|108|) . The outoging solutions are constructed in the usual 
manner from the spherical Bessel and Hankcl functions: 



* = -2iQK rnj e (mnr < )h ( e 1) (milr > )Y; m (Tr/2,0)Y ern (8, l f) 



(114) 



,m even 



or 



* = -2K ]T Y; o (ir/2,O)Y m (0,O) 



+ 4KnJ2 Y ™Y; m (n/2,0)Y im (9,0)j e (mnr < )Im{h£\mnr > )e lm ' p } 



(115) 



m=2,4,6 



By expansion in spherical harmonics, then by the usual Green function construction, the outgoing solution for <J/( n1 ' 
is found to be 



^(nl) = jjinl) + iV (nl) = _ ^m Q Vy ^ 



17- 1 (7r/2,0)r /l _i(fl,0) 



-lif 



21+1 ) 



+n J2 (™ + l)^: m (7r/2,O)^„ l (0,O)^((m + l)r!r < )^ 1) ((m + l)r!r > )e J 

m— odd, 7^—1 



(116) 



' odd 



+ n V (m + l)r/ m (7r/2,0)r £m (6»,0)^((m+l)Or < )[cosm^j £ ((™+l)fir > ) — sinm^n^((m+l)fir>)] (117) 



m=l ,3,5, 



-n ^ (m+l)y^ 2 (7r/2,0)y Vm _ 2 (M)M™+^ 

m— 1,3,5, ■■■ 
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and 



v (ni) = _ 327rmot , 7 ^ 



£ odd 



cos (f r< 



(2£+l)r e + 



T y / % 1 (7r/2, ojy/.-i^.o) 



n 53 (m+l)y^(7r/2,0)y <m (e,0)^((m+l)fir<)[cosm<pn < ((m+l)fir>)+sinm^^((m+l)nr>)] (118) 

m=l,3,5"- 



+fi 53 (m+l)y^_ 2 (7r/2,0)y^_ 2 (M)M™^ 

m=l,3,5- ■ ■ 



Similarly, the outgoing solution for \I/( 22 ) is found to be 

$ ( 2 2) =f/ ( 22 ) +iy (22) = 167rmoW 2 7 J- 



r/_ 2 (7r/2,o)y/,_ 2 (fl,o) 



2^ + 1 ; r £+! 



(m+2)y^ m (V2,O)y, m (0,O)i,((m+2)r!r < )^ 1) ((m+2)r!r > ) 



m— even, 7^ — 2 



(119) 



or 



C/( 22 ) = I6imi v 2 j ^ 



^ even 



— cos 2(p r< 
{2t+l) ri +1 



y/_ 2 (7r/2,o)y,_ 2 (0,o) 



+ n J! (m+2)y/ m (7r/2,O)y£ m (0,O)^((m+2)nr<)[cosm^n < ((m+2)nr>) + sinrrup je((m+2)ilr > )] (120) 



m=0,2,4-. 



+fi 5] (m+2)y^_ m _ 4 (7r/2, 0)y,- m - 4 (6>, 0)^((m+2)Or<)[cos ((m+4)p) n/((m+2)fir>) + sin ((m+4)p) j/((m+2)fir>)] 

m=0,2,4--- 



and 



^ (22) = 167rm w 2 7 53 



^ even 



(|^y^|r^- 2 (V2,o)y ; - 2 ^o) 



- fi 53 (m + 2)y/ m ( 7 r/2,0)y £m (^0). ?£ ((m+2)r!r < )[sin mtp n£((m + 2)f7r>) — cosm</?^((m + 2)£lr > )] (121) 

m=0,2,4"- 



5I( TO + 2 ) y £*-m-i( 7r /2,0)y,^„^(6',0)j £ ((TO+2)fir < )[sin((TO+4)^) n/((m+2)fir>) - cos ((m+4)p) j/((m+2)fir>) 

m=0,2,4--- 



The expressions for U {nl \ V {nl \ U {22 \ and y( 22 ) have been given as sums only over nonnegative values of m through 

P(-x)=Mx)(h? 



the use of the relationship x)h^\— x) = ji(x){h!p (x)) 
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IV. LINEARIZED GRAVITY IN ADAPTED COORDINATES 
A. Field equations in adapted coordinates 

The field equations to be used for computation are the source- free forms of Eqs. 182(1 - l|88|l. The first four of these 
equations are 

□$ (afc) = (122) 
for (ab) — (00), (nO), (20), and (nO). In real form the last three of these equations are 

DU (ab) + 2kQ 2 d v V (ab) + k 2 Q?U (ab) = (123) 

uv {ab) ~ 2kn 2 d v u (ab) + k 2 n 2 v [ab) = o (124) 

with k = 1 for (ab)=(nl), or (21), and k = 2 for (ab) = (22). 

The form of the scalar d'Alembertian in adapted coordinates is given in Eq. (8) of Ref. [l2j, as 

a 2 * d 2 ^ d 2 t> d 2 ^ , a 2 * a 2 * 



= A xx — — + Aee -rr^ + A$$ + 2A xB — — + 2A X $ — — + 2A e * 



xx d x 2 1 do 2 1 "* w a$ 2 1 1 1 aea$ 

a* a* a* 



<9x ae 9$ 

and the adapted-coordinate form of d v in Eq. (27) of Ref. [T^. 

9 ^ e a* , ^ v a* 



B x- + B e^ + B^, (125) 



Here the ^4, -B and T coefficients are known, real functions of Xi 9, $ that are given explicitly in Appendix A of 
Ref. ^2 an d are repeated in the appendix of the present paper. 

B. Inner boundary conditions in adapted coordinates 

To express the inner boundary conditions in terms of the adapted coordinates, we approximate 
r 1 + 7 W = {Z- a) 2 + X 2 +Y 2 + 7 2 v 2 X 2 = 

[l+ 7 Vsin 2 29cos 2 $] ^+0{ X 6 /a 4 ), (127) 
and we write Eqs. I(89|l -194 (l . for the x — > limits of the the fields, as 

*W = 4m 2a72 (128) 

X 2 v 1 + 7 2 w 2 sin 2 26 cos 2 $ 

*W> = 4m - 2072 4= (129) 

X 2 V 1 + l 2v2 sm2 29 cos 2 $ v3 

*™ = 4 TO0 - 2fl ^„ (-*) (130) 



2a 7 2 


xVi- 




26 cos 2 






2a7 2 






xVi- 




29 cos 2 






2«7 2 






xVi- 


|_ ^2^2 gm 2 


29 cos 2 






2a 7 2 






xVi- 


9 9*2 


29 cos 2 






2a7 2 






xVi- 


9 9-2 


29 cos 2 





W nl > = 4m I ^sgn[cos9] (131) 

X V 1 + TV sin 29 cos 2 $ 

t/ (22) = „ , s 2a72 5 f-^) (132) 



^(nO) = jj(21) = y(21) = jj(nl) = y-(22) = q ^33) 
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C. Outer boundary conditions in adapted coordinates 

The outer boundary conditions in adapted coordinates follow from Eqs. (|95|l - l|98|) with r replaced by \ and d v 
replaced by the expressions in Ecl (|126j) . For the fields $(°°), $("°) and *( 20 \ the conditions are just those 

used for the scalar field in Ref. [l^- For the complex fields, the conditions are a modification of Eqs. (|99|l and (|l(J0f) . 
and are the tensorial equivalent of the adapted-coordinate outer boundary condition l|ll|) used for scalar fields; aside 
from corrections of order (x/a) 2 these results are: 

I|(^,). ± ,(^» +r e^ + r .^ +r ^), 

where 

k = 1 for (ab) = (nl), (nO), (21) k = 2 for (ab) = (22) . (136) 



V. NUMERICAL RESULTS 

The field equations in adapted coor dinates for *(°°), *( n °) and *( 20 ) are defined by the homogeneous wave 

equation (|122fl and by the form of the d'Alembertian in Eq. (|125fl with the coefficients in AppendixlX! These equations 
are subject to the inner boundary conditions in Eqs. (|128fl - l|130|) applied at some smallest value x m i n of x- The 
equations must satisfy the outgoing or ingoing outer boundary conditions Eq. (|llf) applied at a largest computational 
value Xmax °f X- I n the linear theory each of the scalar-like fields <I'( rm ), ifK 00 ^ ^r(0°) j \Jj(20) j s completely decoupled 
from every other fields, both in the field equations and in the boundary conditions. Furthermore, the field equations, 
and boundary conditions have precisely the same forms as those for the scalar field problem. The computational 
problem, therefore, is precisely that of Ref. [T^ where it was shown that the computed solution agrees accurately with 
the numerically evaluated series solution in Eq. (|115(l . 

In the computation of the complex fields v]/( nl ) and \I>( 22 ), features arise that are different from those in Ref. [l2T |. 
Here the field equations (112 31) and (J124JK and the outer boundary conditions (iTTI . couple the real and imaginary parts. 
[This would apply also to the complex field \Ia 21 ), but due to the inner boundary condition \I>( 21 ) = in Eq. (|133|l . 
the field \I>( 21 ) must be be identically zero in linearized theory] 

The series solutions for and are given in Eqs. ljTT7|) - (|TT%j) and those for and in Eqs. [jT2U|) 

- I|121|) . Figures 0] and [3] give a comparison, for two different source velocities, between these series solutions and the 
solutions of the eigenspectral method: adapted coordinates and multipole filtering based on the modified multipolcs 
appropriate to the discrete angular operator. 

A feature that stands out in the figures is the disagreement at small x between the series and eigenspectral solutions 
for and for U {22 \ The failing here is in the convergence of the series solutions in Eqs. (|1 18|) and (|120|) . The 

y( nl ) a nd U^ 22 ^ fields diverge at x — > 0, hence these series converge very slowly at small x- Numerical experiments 
summing very large numbers of terms and evalutating them with arbitrary precision arithmetic confirm that the series 
solutions for and [/( 22 ) in Figs. 0] and [5] have large errors. The series for LA" 1 ) and V^- 22 \ on the other hand, 

are convergent at x = and show good agreement with the small- x form of the eigenspectral solutions. 

The series solutions are highly accurate for % > 1, so the differences between the series solutions and the eigenspectral 
solutions for large x are an indication of the limitations of the eigenspectral method. Those differences are larger for 

> 0.4 than for fi > 0.3. Computations (not presented here) for fi > 0.5 show significantly larger error. The origin 
of these errors is the relatively coarse computational grid used, and the limited number of multipoles used in the 
multipole filtering. (This is equivalent in the computation to the coarseness of the angular grid.) The error is simply 
that due to truncation error, and is expected. As Q is increased the gradients of the fields increase and truncation 
error becomes more important. These limitations are imposed by the fact that the computation was done on a 2GB 
RAM workstation. Greater accuracy, and hence higher velocity, would be possible on larger machines, 

To test whether our numerical techniques might be sufficient for the next steps in our program, we introduced a 
simple nonlinearity into the equations of Sec. IIVI 



□F = D$W = k—„ (137) 

H 2 + a 2 S 
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FIG. 4: A comparison, for a£l = 0.3, of outgoing linearized gravitational fields computed by series summation (solid curve), 
and by the solution of the coupled partial differential equations of the "eigenspectral" method of the PSW approach (dashed 
curve). The fields are shown along a line (the O = line) outward through the source. For these computations, the grid in \, 
O, $ was 1500, 16, 32, respectively. The entire angular space, < < 7r, < $ < 2n, was used. The inner boundary was at 
Xmin = 0.1a and the outgoing boundary condition was imposed at Xmax = 30a. The multipole filtering included all modes up 
through the octupole, approximately I = 3. (The actual discrete eigenvalues differ slightly from integer values.) See the text 
for a discussion of the numerical limitations of the series summation. 
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FIG. 5: A comparison, for afl = 0.4, of outgoing linearized gravitational fields computed by series summation (solid curve), 
and by the solution of the coupled partial differential equations of the "eigenspectral" method of the PSW approach (dashed 
curve). The details of the computations are the same as those for Fig. |1| 



with S defined by 



b = f] h, \ njrr) a p 



where h na = (n • h) a . 

With the notation of Eqs. l|54fl this can be written as 



(138) 



(139) 



since $W and are "helical scalars," i.e. , functions only of corotating coordinates, for (ab) = (nn) or (nO), we 

have 

^^(06)^(06) T ,Aa^(o6)^(ab) _ 
'/ ,A ,a I ,\ ,u 



n 2y(ab)y(ab) + . y(ab)y(ab) + (yQ . ye) + (y$ _ y$) $ 



r ,( ab ),T,( ab ) 



(140) 
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in which 







ve 


ve 


v$ 


v$ 



Q_ 

x 2 
Q_ 

x 4 



= 2 



Q + a 2 + x 2 cos(2 9) 
X 4 sin 2 (2 6) 



Q = v / a 4 + 2o 2 x 2 cos(26) + x 4 
For we have (\£r( nl ) = e - int ^ii nl ) which results in 



(141) 
(142) 

(143) 
(144) 



- or 



$(nl) ^C" 1 )^* _ ipC™ 1 ) Aj/O^A* 

(Vx • vx) (% x) ) * + (ve ■ ve) ($ g 1} V + ( v$ ■ v $) * 



^ ( ^ 4) ^ ( ™ 4) ^ 



{/(nl) 



2ft 2 



(Vx • Vx) 



(e') 2 +(^ 1: 



(ve • ve) 



«)) 2 + « 1} ) ; 



ft 2 



(V* • V$) 



+ (^ nl) ) 



(145) 



In the modified "theory" represented by Eq. (jl37[l the equations that determine $( n0 ) and remain unchanged, 

so these fields are found with linear equations. The nonlinear occurrence of vpC™™) i n S, however, means that vIm™™) 
solves a nonlinear equation. 

The solution to this nonlinear problem requires the iterative techniques that have previously been used for nonlinear 
scalar models in Ref. ^| • Results from the applications of these iterative methods to the modified theory are displayed 
in Fig. El for various values of the nonlinearity parameter k. For all values of k the inner boundary conditions are 
taken to be those of the linear problem and are imposed at Xmin =0.1; for all models outgoing boundary conditions 
are imposed at Xm&x = 30. A comparison, in that figure, with the k = linear solution demonstrates that the cross- 
coupling of fields and nonlinearity has a strong effect on the solution to the toy model, changing the amplitude of the 
vj/(™™) waves by an order of magnitude. The purpose of this nonlinear toy computation is not to extract physics, but 
simply to suggest that the iteration techniques previously developed will be adequate at least for a range of nonlinear 
models. 



VI. SUMMARY AND DISCUSSION 

This paper has laid out the infrastructure for describing helically symmetric linearized gravity, and more generally 
for describing helically symmetric second-rank tensors in a flat background. The fields have been written in terms 
of "helical scalars," i.e. , functions only of corotating coordinates (equivalently, labels on the curves of the Killing 
congruence). The paper has, furthermore, shown how to formulate computational problems in linearized general 
relativity in terms of these helical scalars. The field equations, sources, and inner and outer boundary conditions have 
been written in this format in Sec. IIIII 

A welcome feature of the formulation in Sec. IIIII is the extent of the separation of "components" of the field. In 
the starting point, the Lorentz-gauge field equations = —IGttT^, a separate equation is satisfied by each 

component h^ v with respect to an inertial Minkowski basis. That attractive feature cannot be taken over directly 
to the helically symmetric problem since the components with respect to the inertial basis are not helical scalars. It 
turns out, however, that the formulation of the linear problem in helical scalars leads to four real fields vj>(™°) ; 
^r(oo) anc j vjjr(20) ; anc j three complex fields \p( 21 ) and \p( 22 ), and that these four real and three complex fields 

are not mixed by the field equations or the inner or outer boundary conditions. The only mixing in the problem is 
between the real and imaginary parts of the complex fields. 

Since the three corotating coordinates are general, this paper has also presented the explicit formulation of the 
computational problem in the adapted coordinates that were found in Ref. [12j to be extremely useful. The numerical 
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FIG. 6: The field ^ nn '(x) f° r the toy nonlinear model described in the text. Curves are marked with the value of the 
nonlinearity parameter k; the parameter H was set to unity and the models were run on a 1500 x 16 x 32 grid for the full 
X, 6, $ space (i.e. , no symmetries were used), with Xmin/a = 0.1, and Xmax/a = 30. In the multipole filtering, modes through 
octupole were included for all fields. Comparison with the k = linear solution shows the importance of nonlinear effects. 



results in Sec. of the present paper demonstrate that the numerical challenges presented by linearized gravity are 
the same as those for the linear scalar model in Ref. 12}]. Indeed, the numerical problem for the four real fields is 
exactly the same as that for the linear scalar field. For the complex fields, the new features are minor modifications of 
the boundary conditions and mixing of the real part and of the imaginary part of each of the complex fields. These 
new features do not appear to present any new numerical difficulty, and in fact no difficulty was found. Furthermore, 
a trial with a toy nonlinearity suggests that there are also no new problems in dealing with nonlinear terms, except 
those of complexity. 

Although the present paper deals almost exclusively with linearized gravity, the infrastructure developed here is 
more widely applicable. Our next step in the PSW program is to solve the post-Minkowskian problem for the orbiting 
point masses. The equations to be solved in this method have the same operator /i ^ a on the left hand side, but 
have "sources" quadratic in h^ u ^ a on the right. (There is also a second derivative on the right multiplied by an 
undifferentiated /i M „; this term can be treated, like the others, as a perturbation, or it can be moved to the left to 
modify the principal part.) In the usual spirit of a post-Minkowski approximation, we could solve first for the first- 
order fields and treat them as known sources. Alternatively, in a numerical approach, we could treat the equations 
as a given nonlinear problem. Either way, the formalism developed in the present paper goes over directly to the 
post-Minkowskian problem. Again, the fields can be described with four real and three complex helical scalars, and 
the equations and boundary conditions follow from the simple relations between the helical scalars and their inertial 
equivalents, i.e. , relations like those in Eqs. Ij64[l - l|7U|) . 

The full Einstein equations can also be viewed as a higher order extension of the post-Minkowski equations. In 
principle, the only change from the post-Minkowski problem is the inclusion of terms of all orders on the nonlinear 
right-hand side of the field equations. Again the formalism developed in the present paper should suffice for the 
description of the problem, and should be convenient. In particular, helical scalars of the background Minkowski 
space will be helical scalars of the full metric. In practice, new problems will arise. One is the question of the 
relativistic Kepler's law: what is the appropriate relationship relating the source strength (encoded in inner boundary 
conditions), the coordinate separation of the sources, and the parameter f2? The post-Minkowskian approximation, 
for which the answer is known, will help clarify how this is to be handled in full general relativity. A more subtle 
question is whether it is justified to use a formalism based on weak-field structures to describe strong gravitational 
fields. For highly curved spacetimes do Minkowksi-like coordinates exist with which we can use the formalism of the 
present paper? Possibly relevant to this is the fact that in our computations we can impose inner boundary condistions 
at some distance from the sources, so that the effects of strong fields can be somewhat controlled. 

In any case, the relative simplicity of the description presented here, along with the absence of any new computa- 
tional difficulties (that is, difficulties not present in the scalar problem) is a reason for optimism that the next steps 
can be taken reasonably quickly. 
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APPENDIX A: COEFFICIENTS FOR ADAPTED COORDINATES 

The adapted-coordinate coefficient A, B, and T are listed here. Derivations are provided in Ref. [l2|. 
As in Ref. 12) the coefficients needed for the wave operator in Eq. I|125|) are written in the form: 

A xx = — - Sl 2 A xx (Al) 

A 

Aee = ^-n 2 A ee (A2) 

X 4 sin 2 (2 6) V ' 

A xe = -n 2 A xe (A4) 

A e * = ~VL 2 A e<s , (A6) 
B x = a ^±p- - n 2 B x (A7) 

B B = V Q^ 2 +X 2 cos(2B) (Q-a*l _ ^ 

x/Q-a 2 - x 2 cos(2 9) X 4 
B$ = -n 2 B^, (A9) 

where Q is given in Eq. (|144|) . 

The expressions multiplied by fl 2 are: 

- a 4 sin 2 (29) cos 2 $ tk 
A xx = W ( A1 °) 

A 

cos 2 $ Ty 2 + a 2 cos(29)l 2 
Aee l - 4 (All) 

t . , Q + a 2 + x 2 cos(29) .. . 

^ =Sm %- a 2 -X 2 cos(29) (M2) 

- a 2 [x 2 + a 2 cos(29)] sin(29) cos 2 $ 

A x0 = -3 (A13) 

A 

a 2 [Q + a 2 + x 2 cos(29)] sin $ cos $ 
A x<s> = -3 ( Al4 ) 

A 

sin($) cos($) [a 2 + x 2 cos(2 9) + Q] W 2 + a 2 cos(2 9)1 

Ae* = 4 . t^tv; ( Al5 ) 

X sm(2 9) 

a 2 [cos 2 ($) {3a 2 cos 2 (29) - Q - 2a 2 + x 2 cos(29)} + Q + a 2 + x 2 cos(29)] 

B x = -3 (A16) 

A 
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(3Q + a 2 + x 2 cos 26) sin($) cos($) 
S * - Q-a2_ x 2 cos2e ( A17 ) 

xVQ-a 2 -X 2 cos(26) 

where 

c = a 2 x 4 cos(2 8) + 2 a 4 X 2 + 4 a 6 cos(2 6) + 4 a 4 x 2 (cos(2 8)) 2 - 4 a 4 Q cos(2 9) - 2 a 2 Q X 2 - X 6 (A19) 

d = x 4 (a 2 cos(2 6) + x 2 ) ■ (A20) 
The coefficients needed in Eq. I|126|) a to express the Sommerfeld boundary condition in adapted coordinates are: 

r x = (z% - X d X\ = a2cOS$sin(26) (A21) 
V dX dzj X 

r° = (z™ - x°2) = -^(« 2 ^(2e) + x 2 ) 

V dX dZj X 2 

r , = /-M _ = _ X 2 sin^in(29) 

V aZ/ -a 2 -x 2 cos26 + Q 

Note that there were errors in the expressions given for the Ts in Ref. [l2|. 
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